Technical  Research  Report 

Formation  Dynamics  under  a  Class  of  Control  Laws 


by  Fumin  Zhang,  P.S.  Krishnaprasad 


CDCSS  TR  2002-4 
(ISR  TR  2002-25) 


r 


CENTER  FOR  DYNAMICS 
AND  CONTROL  OF 
SMART  STRUCTURES 


The  Center  for  Dynamics  and  Control  of  Smart  Structures  (CDCSS)  is  a  joint  Harvard  University,  Boston  University,  University  of  Maryland  center, 
supported  by  the  Army  Research  Office  under  the  ODDR&E  MURI97  Program  Grant  No.  DAAG55-97-1-0114  (through  Harvard  University).  This 
document  is  a  technical  report  in  the  CDCSS  series  originating  at  the  University  of  Maryland. 


Web  site  http://www.isr.umd.edu/CDCSS/cdcss.html 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2002 

2.  REPORT  TYPE 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

Formation  Dynamics  under  a  Class  of  Control  Laws 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Army  Research  Office, PO  Box  12211, Research  Triangle  Park, NC, 27709 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

The  original  document  contains  color  images. 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF  PAGES 

9 

19a.  NAME  OF 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Formation  Dynamics  Under  A  Class  of  Control  Laws1 
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Abstract 

A  system  of  two  earth  satellites  is  analyzed  as  a  con¬ 
trolled  mechanical  system.  The  orbit  of  an  earth  satel¬ 
lite  can  be  represented  by  a  point  in  the  vector  space  of 
ordered  pairs  of  angular  momentum  and  Laplace  vec¬ 
tors.  Control  laws  are  obtained  by  introducing  a  Lya¬ 
punov  function  on  this  space.  Formations  of  two  satel¬ 
lites  are  achieved  asymptotically  by  the  controlled  dy¬ 
namics. 


1  Introduction 

There  is  increasing  interest  in  the  behaviors  of  a  cluster 
of  relatively  small  satellites  in  space.  Such  satellites  are 
usually  inter-connected  by  wireless  radio  or  laser  links 
for  communication.  By  keeping  a  cluster  of  such  satel¬ 
lites  in  a  certain  geometrical  form,  one  can  acquire  ben¬ 
efits  for  scientific  observations.  The  information  shar¬ 
ing  across  the  cluster  will  allow  the  satellites  to  work 
cooperatively  to  perform  tasks  impossible  or  difficult 
for  a  single  satellite.  Compared  to  a  single  satellite 
providing  the  functionality  of  a  cluster,  a  member  of 
a  cluster  can  be  smaller/lighter.  Building  and  launch 
costs  will  then  be  reduced.  In  addition,  A  cluster  can 
also  be  reconfigured  according  to  different  mission  goals 
or  in  the  case  that  a  member  of  the  cluster  fails. 

The  size  and  shape  of  a  cluster  or  formation  are  usually 
determined  by  the  required  functionality  of  the  forma¬ 
tion.  With  J2  and  higher  order  terms  of  the  earth’s 
gravitational  field  ignored,  the  solution  of  the  Kepler 
problem  tells  us  that  a  satellite  will  track  an  elliptical  or 
circular  orbit.  To  determine  the  orbit  of  each  satellite 
in  a  cluster  such  that  the  size  and  shape  of  the  cluster 
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is  kept  unchanged  over  sufficiently  long  period  is  not  a 
trivial  problem.  Because  the  amount  of  fuel  on  board 
is  limited  for  each  satellite,  one  needs  to  find  a  set  of 
orbits  that  demand  the  minimum  control  effort.  This 
problem  becomes  more  challenging  when  the  effect  of 
J2  and  other  disturbances  are  considered.  In  [3]  the 
authors  proposed  a  ring  of  evenly  distributed  satellites 
on  the  same  circular  orbit  for  communication  purpose. 
The  stability  of  such  ring  is  proved  in  [4].  In  [2]  and  [5] 
the  investigations  of  Clohessy- Wiltshire  equations  re¬ 
vealed  possible  formations  with  constant  apparent  dis¬ 
tributions.  The  effects  of  perturbations  are  calculated 
and  possible  station  keeping  strategy  are  proposed.  In 
[7],  the  authors  proposed  that  in  the  presence  of  J2,  a 
set  of  constraints  on  the  orbital  elements  shall  be  satis¬ 
fied  to  prevent  the  orbits  from  drifting  apart.  However, 
extra  station  keeping  is  still  necessary  due  to  the  com¬ 
plicated  nature  of  the  disturbances.  The  adjustment 
can  be  performed  periodically  when  the  drifting  error 
exceeds  certain  threshold. 

The  initialization  of  a  formation  is  another  important 
problem.  The  whole  cluster  can  be  launched  together 
by  a  space  shuttle  or  rocket.  Satellites  will  be  first 
placed  in  a  parking  orbit  before  transferring  to  the  final 
orbits.  The  final  orbits  must  be  such  that  the  forma¬ 
tion  will  be  achieved.  The  orbital  transfer  can  happen 
individually.  One  idea  is  to  develop  a  Lyapunov  func¬ 
tion  which  will  achieve  its  minimum  when  correct  orbit 
is  reached.  In  [6],  a  Lyapunov  function  is  expressed  as  a 
quadratic  function  of  the  differences  of  orbital  elements 
between  current  orbit  and  the  destination  orbit.  In  [1], 
the  authors  proved  that  an  elliptic  orbit  can  be  repre¬ 
sented  uniquely  as  a  point  on  the  linear  space  formed 
by  the  angular  momentum  vectors  and  Laplace  vectors. 
A  Lyapunov  function  was  built  naturally  from  the  Eu¬ 
clidean  metric  on  this  space.  It  has  been  suggested  that 
compared  to  Hohmann  transfer,  this  approach  will  con¬ 
sume  less  fuel. 

In  this  paper,  we  will  show  that  cooperative  orbit  trans¬ 
fer  is  possible  in  the  setting  of  pairs  of  satellites.  The 
satellites  in  a  cluster  can  make  the  transfer  together 
and  the  cluster  relationship  can  be  established  asymp¬ 
totically.  The  same  technique  can  be  used  to  perform 


station  keeping  task  when  the  relative  position  errors 
of  the  satellites  exceed  a  certain  threshold.  By  extend¬ 
ing  the  method  in  [1]  to  multi-satellite  case,  we  will  de¬ 
sign  a  Lyapunov  function  and  show  that  correct  cluster 
relationship  between  orbits  can  be  established  by  the 
controlled  dynamics  which  will  minimize  this  function. 

In  section  (2),  we  will  review  the  Kepler  two  body  prob¬ 
lem  and  introduce  the  concept  of  the  shape  space  of 
elliptic  Keplerian  orbits.  In  section  (3),  the  Lyapunov 
function  for  a  pair  of  satellites  is  developed  and  the  case 
of  two  satellites  making  a  cooperative  orbit  transfer  is 
studied.  Simulation  results  are  given  in  section  (4). 


2  The  Kepler  two  body  problem 

For  a  system  of  two  small  satellites,  one  can  make  the 
following  assumptions:  (a)  The  gravitational  attraction 
between  the  satellites  can  be  omitted,  (b)  The  total 
mass  of  the  satellites  satisfies  m i  +  m2  M  where  M 
is  the  mass  of  the  Earth.  Under  these  assumptions,  the 
three  body  system  can  be  approximated  by  two  uncou¬ 
pled  two  body  problems.  Each  of  the  two  body  prob¬ 
lems  can  be  further  simplified  to  a  one  center  problem 
with  the  center  of  Earth  being  the  center  of  mass. 

The  gravitational  potential  function  V  depends  only  on 
the  distance  ||</||  of  the  satellite  from  the  center.  Then 
we  have 

mq  =  — V  V  (1) 


Let  p  =  mq  be  the  momentum  of  the  satellite.  The 
angular  momentum  l  =  q  x  p  and  the  energy  W  = 
im||4||2  +  V  are  integrals  of  motion  (i.e.  they  are  con¬ 
served).  Since 


and 


=  q  X  p  +  q  X  p 
=  0  +  q  x  mq 

q  dV 

=  —q  X  m-rr- 


=  0 


W  =  q  ■  mq  +  V 

=  -q-\7V  +  V  =  0 


(2) 


(3) 


If  we  make  further  assumptions  that  the  shape  of 
the  earth  is  a  perfect  homogeneous  ball,  then  V  = 
— m/i/||g||  where  p  =  kM  and  k  is  the  gravity  con¬ 
stant,  the  Laplace  vector  A  =  p  x  l  —  rri2pj^  is  also 
conserved  given  q(t)  ^  0  for  all  t  >  0.  To  see  this, 


4  -7  2  d  r  1  \ 

.4  =  pxi-,n„-(M) 


2..,  q  {q-q)q 


=  mq  x  (q  x  p)  —  m  m(fT  “ 


)  (4) 


Notice  that 


Thus 


mq  = 


— mp 


q 


IMI3 


mq  x  (q  x  p)  =  — mp 


(■ q-p)q 

IMI3 


(5) 

(6) 


Comparing  (6)  with  (4)  we  have  A  =  0 

Knowing  l  ,  W  and  A,  we  have  seven  integrals  for  the 
two-body  problem.  They  are  not  all  independent  be¬ 
cause  there  are  two  relations  connecting  them.  They 
are: 


A  - 1  =  0 

||x4||2  =  m4 p2  +  2mW  \\  l  \\2  (7) 


The  space  of  ordered  pairs  (l,  A)  is  R3  x  R3  on  which 
we  define  the  metric: 

ddluAJ,  0 h,A2 ))  =  (||  h  -  l2  f  +  Mr  -  A2  ||2)^ 

(8) 

Let  P  denote  the  phase  space  of  the  satellite,  define  a 
mapping  7r  :  P  — »  R3  x  R3,  ( q,p )  i— >  (l,  A).  Let  the  set 
Ee  be  defined  as 

Xe  =  {(q,p)eP\W(q,p)<  0,1^0}  (9) 

and  let  the  set  D  be  defined  as 

D  =  {(l,  A)  e  R3  x  R3  I  A  ■  l  =  0,  M  0,  \\A  ||  <  m2p} 

(10) 

In  [1],  the  authors  proved  the  following  results: 

Theorem  2.1  (Chang-Chichka-Marsden)The  follow¬ 
ing  hold: 

1.  Ee  is  the  union  of  all  elliptic  Keplerian  orbits. 

2.  7r(Ee)  =  D  and  S e  =  7r_1(U). 

3.  The  fiber  7r~1(l,A)  is  a  unique  (oriented)  elliptic 
Keplerian  orbit  for  each  (l,  A)  £  D. 


The  mapping  7r  is  a  continuous  mapping  because  l  and 
A  are  continuous  with  respect  to  ( q,p ).  Furthermore, 
the  following  corollary  hold. 

Corollary  2.2  7r-1(/r)  is  compact  for  any  compact  set 
K  C  D 

Proof:  Since  D  is  a  metric  space,  K  is  compact  implies 
that  K  is  closed  and  bounded.  By  the  continuity  of  7 r, 
7t_1(A")  is  closed.  The  only  thing  left  to  prove  is  the 
boundedness  of  7t_1(A"). 
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K  is  closed  and  bounded  implies  that  there  exist  r o  >  0, 
ri  >  0  and  0  <  r2  <  m2p  s.t. 

?’o  <  IM  II  —  ri  II  A II  <  r2  (11) 

This  is  true  since  we  already  know  ||  A  |  <  m2/i  and 
||  Z  ||  >  0.  Because  ||  •  ||  is  continuous  on  the  compact  set 
AT,  we  declare  that  ro  and  r2  exist  and  can  be  achieved. 
The  existence  of  r\  is  based  on  the  fact  that  ||  l 1|2  + 
||  A  || 2  is  bounded. 

Without  loss  of  generality,  assume  the  satellite  has  unit 
mass.  From 


IM II2  =  f*2  +  2ir  ||  ( ||2 

(12) 

We  get 

(13) 

2  II Z  f 

Hence, 

0  T2  -  rl  <  ijy|  <  P2 

2r?  -1  '-2r2 

(14) 

Let  a  denote  the  semi-major  axis  of  an  elliptic  Kep- 
lerian  orbit  and  e  denote  the  eccentricity.  It  is  well 
known  that 


(15) 


Thus  a,e  are  bounded.  For  any  ( q(t),p(t ))  G  7t_1(A"), 
on  an  elliptic  orbit  q(t)  is  bounded  below  by  a(l  —  e) 
and  bounded  above  by  a(l  +  e).  Furthermore,  because 
||  1 1|  =  ||  p(t)  ||  ||  q{t)  ||  sin  0,  ||  p(t)  ||  is  bounded  given  the 
fact  that  ||  Z  ||  7^  0.  Hence  ( q(t),p(t ))  is  bounded  which 
implies  that  tt~1{K)  is  bounded.  ■ 


Proof:  Use  the  definitions  of  (l,  A,  Ee,  D,n,  d) 

and  apply  Theorem  (2.1)  and  corollary  (2.2)  to 
(li,  Ai,Yjei,  Di,  ni,  di)  for  i  =  1,2.  The  statements  are 
immediately  proved.  ■ 


3  Control  Strategies  to  Achieve  Formations 

Based  on  the  distance  function  on  the  shape  space  of 
elliptic  Keplerian  orbits,  a  Lyapunov  function  can  be 
introduced  on  the  phase  space  for  the  Hamiltonian  sys¬ 
tem  of  two  satellites.  By  applying  LaSalle’s  invariance 
principle,  we  prove  that  the  system  can  be  driven  to  an 
invariant  manifold  of  the  phase  space  where  formation 
is  achieved. 

The  system  of  two  satellites  is  a  Hamiltonian  (Control) 
system  with  H  =  y~|;2  ZZ,  where 

Hi  =  ^-\\Pi\\2-V(\\qi\\)  (17) 

Zrrii 

For  each  satellite  the  dynamics  is: 

dH,  q, 

— ~ - 1 -Ui  =  -m,ipT — yr  +  Ui 

oqi  lift  || 3 

^  —  (18) 

dpi  rrii 

where  i  =  1,2.  Here  u,  are  controls. 

Using  the  notations  at  the  end  of  the  last  section,  let 
Ji{diP)  =  {h.M),  Mq,p)  =  (h,A2)  and  J(g,p)  = 
(Ji,  J2).  Define  a  Lyapunov  function  as 


Pi  = 

ft  = 


In  the  case  of  two  satellites,  let  (k,  A-i ,  Pi ,  Eei,  Di:  Tri:  di) 
denote  the  corresponding  objects  defined  for  the  ith 
satellite.  Let  P  =  Pi  x  P2  ,  q(t)  =  (qi(t),  q2(t)),  p(t))  = 
(pi(t),p2(t)).  Let  Ee  =  Eei  x  Se2-  Let  D  =  D\  x  D2 
and  l  =  (li,l2),  A  =  (Ai,  A2).  Let 

d{(Ji,Ai),  (f2,  A2))  =  d\  +  d\  (16) 

and  7r  =  7Ti  x  7T2.  We  define  the  shape  space  of  ellip¬ 
tic  Keplerian  orbits  to  be  the  set  D  with  the  distance 
function  d. 

Proposition  2.3  The  following  hold: 

1.  Ee  is  the  union  of  all  pairs  of  elliptic  Keplerian 
orbits. 

2.  7r(Ee)  =  D  and  Ee  =  7f_1(U). 

3.  The  fiber  7 r_1(Z,  A)  is  a  unique  pair  of  (oriented) 
elliptic  Keplerian  orbit  for  each  (l.  A)  €  D. 

f.  7r~ 1  (AT)  is  compact  for  any  compact  set  K  C  D 


V(J(q,p))  =  ^[||Zi  —  l2  —  5l\\2  +  || Aa  -  A2  —  <5A||2 
+  ||  Zi  —  Id  ||2  +  ||  An  —  Ad  ||2]  (19) 

Here,  51  ,  5A  are  constant  vectors  which  specify  the 
desired  final  difference  between  (li,  A\)  and  ( l2,A2 ).  Id 
and  Ad  are  also  constant  vectors  which  specify  the  ref¬ 
erence  orbit  of  the  formation.  It  can  be  verified  that  V 
is  continuously  differentiable  and  bounded  from  below 
by  0. 


To  investigate  the  derivation  of  the  control  law,  we  can 
rewrite  (18)  as: 


x  dH 

P  =  +  u  =  B[q,p)  +  u 


Q  = 


dH 

dp 


(20) 


and  J(ftp)  is  conserved  when  u  =  0.  Then  we  have 
dV  dJ  x  dV  dJ  ^ 

dV  dJ  x  dV  dJ _ dV  dJ 

dJ  '  dqq+  dJ  '  dpB^q,P^  +  dJ  '  dpU 


By  letting 


dV  •  dV  dJ 

-  dJ  ■  J\u=  0  +  dJ  ■  g~U 

dV  dJ  d J  TdV 

~dl'~d§u  =  ^  ~dl'  u 

(21) 

By  letting 

dJTdV 

“  =  aJ> 

(22) 

i.e 

,  ^,dV  dJ 

(23) 

where  A  =  diag(Xi)  and  A;  >  0  for  i  =  1,2,  we  have 
V  <  0  along  the  integral  curves  of  the  system.  In  [1]  the 
authors  used  this  observation  for  the  setting  of  a  sin¬ 
gle  satellite  to  derive  a  Lyapunov  based  orbit  transfer. 
The  following  theorem  based  on  LaSalle’s  invariance 
principle  holds. 

Theorem  3.1  Given  the  system  as  defined  in  (20), 
Let  V  :  P  — >  R  be  a  continuously  differentiable  func¬ 
tion  which  is  bounded  from  below.  Suppose  there  exists 
a  compact  set  LI  C  Ee  on  which  V  <  0.  Let  E  C  LI  be 
the  set  on  which  V  =  0.  Let  M  be  the  largest  invari¬ 
ant  set  in  E.  Then  by  applying  the  control  u  as  defined 
(22),  every  solution  of  the  system  starting  within  Ll  ap¬ 
proaches  M  as  t  — »  oo.  Furthermore,  Let  T  denote  the 
subset  of  Ll  where  u(t)  =  0  holds  for  all  t,  then  M  =  T 

Proof:  It  is  easy  to  verify  that  the  conditions  for 

LaSalle’s  invariance  principle  are  all  satisfied.  Hence 
(q(t),p(t))  — >  M  as  t  — *  oo. 

The  set  T  where  u(t)  =  0  is  an  invariant  set  contained 
in  E.  This  is  because  we  must  have  V  =  u  ■  u  =  0  on 
r.  Hence  T  C  M.  On  the  other  hand,  M  is  invariant 
implies  V  =  0  on  M.  Furthermore,  u(t )  =  0  on  M. 
Hence  M  C  T.  We  conclude  that  M  =  I  .  ■ 

From  this  theorem  we  know  that  the  largest  invariant 
set  is  M  =  finC  where  C  =  {( q,p )  G  Ee|  u(t)  = 
0}.  Then  two  questions  should  be  answered.  The  first 
question  is  whether  a  suitable  compact  set  Ll  exists  in 
Ee.  The  second  question  is  to  determine  C. 

The  following  proposition  gives  the  answer  to  the  first 
question. 

Proposition  3.2  Given 

II  Ml  7^  0 

II  Ad  ||  <  m\g 

Ad -ld  =  0  (24) 

Let  c,Ci,C2,C3,C4  be  positive  numbers  which  satisfy 

ci  = 


C2  = 

\wuw2 

C3  = 

\(m\h-  \\Ad  ||)2 

c4  = 

l-(mlg-\\Ad-8A\\Y 

c  < 

min{c\,  C2,  C3,  C4} 

(25) 

Let  the  set 

G  =  {(T,A)  g  1Z6  x  1Z6\ V(l,A)  <  c }  (26) 

where  V  is  defined  in  equation  (19). 

Then  the  set  Ll  =  n~1(G  Of))  is  a  compact  subset  of 
He- 

Proof:  According  to  proposition  (2.3),  all  we  need  to 
show  is  that  the  set  G  D  D  is  a  compact  subset  of  the 
set  D.  According  to  our  definitions, 

D  =  {(l,  A)  G  R6  x  R%  ■  Ai  =  0,  k  0, 

\\Ai  ||  <  mfp,  i  =  1,2}  (27) 

Obviously,  this  set  D  is  not  a  closed  subset  of  1Z6  xlZ6. 
If  we  let 

K  =  {(?,  A)  GR6  x  R6\k  ■  Ai  =  0,  i  =  1, 2}  (28) 

This  set  K  is  a  closed  subset  of  1Z6  x  1Z6.  The  set  D 
is  a  subset  of  K.  It  is  not  difficult  to  see  that  the  set 
( K  —  D)  is  a  closed  subset  of  1Ze  x  1Z6. 

In  the  next  step  we  want  to  show  that  if  the  value  of  c 
is  chosen  as  proposed,  we  must  have 

G  0  (K  -  D)  =  0  (29) 

Suppose  this  is  not  true.  Notice  that  since  GC\{K  —  D) 
is  a  compact  subset  of  1Z6  x  1Z6,  the  function  V  has  a 
minimum  and  a  maximum  value  on  this  set.  The  maxi¬ 
mum  value  shall  not  exceed  c.  The  minimum  value  can 
be  found  by  solving  a  constrained  minimization  prob¬ 
lem.  One  should  compare  the  unconstrained  minimum 
value  with  the  minimum  value  achieved  on  the  bound¬ 
ary  of  (K  —  D) .  In  this  case,  the  unconstrained  mini¬ 
mum  value  of  V  is  0  which  can  not  be  achieved  since 
(Id,  Ad)  G  TA.On  the  other  hand,  Ci  can  be  achieved 
on  the  boundary  of  (K  —  D )  where  1 1  =  0.  ci  can 
be  achieved  on  the  boundary  where  I2  =  0.  C3  can  be 
achieved  on  the  boundary  where  A\  =  m\p.  C4  can  be 
achieved  on  the  boundary  where  A 2  =  m^g.  However, 
if  c  <  min{c\,  C2,  C3,  C4},  then  the  maximum  is  less  than 
the  minimum.  We  have  a  contradiction.  Thus  the  set 
G  must  have  no  intersection  with  the  set  (K  —  D) . 

On  the  other  hand,  we  know  that  G  0  K  ^  0.  This 
means  G  D  K  C  D.  The  compactness  of  G  0  K  comes 
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from  the  fact  that  G  is  compact  and  I\  is  closed  in 
U6  x  n6 


Next  we  want  to  characterize  the  invariant  set.  By 
equation  (23),  without  loss  of  generality,  we  let  Ai  = 
A2  =  1,  the  controls  can  be  calculated  as 

dV  dh  dV  dAi 

Ul  dl\  dpi  dAi  dpi 

=  — [(ii  —  h  —  SI  +  h  —  ld)xQi 
Ah  x  ( Ai  —  A2  —  5  A  A  A\  —  Ad) 

+  ((Ai  —  A2  —  5A  A  A\  —  Aci)xpi)xgi](30) 

dV  dh  dV  dA2 
dh  dp2  +  dA2  dp2 

=  (h  —  I2  —  6l)xq2  A  ^2x(Ai  —  A2  —  SA) 
T((Ai  A2  SA)xp2)xq2  (31) 

We  need  the  following  lemma  to  solve  the  equations 
obtained  by  letting  ui{t)  =  0  and  U2(t)  =  0. 

Lemma  3.3  Suppose  a  single  satellite  has  external 
control  u[t)  =  0.  Let  x,  y  be  time  invariant  unknown 
vectors.  Suppose  ( q{t),p{t ))  €  S e,  the  solution  of  equa- 


tion 

xxq  +  Ixy  +  {yxp)xq=  0 

(32) 

is 

H 

II 

s 

II 

(33) 

For  some  a  €  R 


Proof:  Take  the  inner  product  with  q(t)  on  both  sides 
of  equation  (32)  we  get: 

(lxy)-q(t)=  0  (34) 

Since  the  orbit  of  the  satellite  is  elliptic,  q(t)  ^  0.  Fur¬ 
thermore,  q{t)  would  stay  in  the  orbital  plane  lA  which 
is  perpendicular  to  l.  Notice  because  Ixy  is  time  in¬ 
variant  and  q{t)  is  not  0,  we  conclude  that  Ixy  must  be 
either  parallel  to  l  or  vanish.  However,  since  the  rela¬ 
tion  (l  x  y)  ■  l  =  0  must  be  satisfied,  the  only  possibility 
is  to  have  Ixy  =  0.  Thus,  there  exists  a  number  a  £  R 
such  that  y  =  al.  Equation  (32)  can  be  simplified  to 


(a:  +  alxp)xq  =  0 

(35) 

Since 

pxl  =  A  +  m2p 

Ik  II 

(36) 

we  now  have, 

{x  —  aA)xq{t)  =  0 

(37) 

Again,  since  q{t)  is  tracking  an  elliptic  orbit,  from  the 
fact  that  x  and  A  are  time  invariant,  we  must  have 
x  =  aA.  m 


Problem  3.4  Suppose  two  satellites  are  controlled  by 
(30)  and  (31)  respectively.  If  we  select  Id  ^  0  and 
II  ^  ||  <  mfp,  we  want  to  calculate  the  maximal  set 
C  C  Ee  where  Ui{t)  =  0  and  U2 (f)  =  0  are  satisfied. 

Solution:  From  ui{t)  =  0,  we  have 

()i  —  I2  —  SI  A  h  —  ld)xqi  +  Ax  {Ai  —  A2  —  +  Ai 

— Ad )  +  ((Ai  —  A2  —  8  A  +  Ai  —  Ad)xpi)xqi  =  0(38) 

From  U2{t)  =  0,  we  have 

(l  1  —  I2  ~~  Sl)xq2  +  l2x{Ai  —  A2  —  SA) 

+  ((Ai  —  A2  —  5A)xp2)xq2  =  0  (39) 

We  want  to  solve  these  two  equations  for  h,l2,Ai  and 
A2.  Notice  that  these  unknowns  are  time  invariant. 

By  applying  the  results  in  the  lemma  to  equation  (38), 
we  get 

1 1  —  I2  —  SI  +  h  —  Id  —  ocAi 

Ai  —  A2  —  SA  -(-  Ai  —  Ad  =  cxli  (40) 

From  equation  (38)  we  get 

Ai  —  A2  —  SA  =  (3l2 

li  ~  I2  ~  SI  =  (3A2  (41) 

Hence, 

h  —  Id  =  olAi  —  (3A2 
Ai  —  Ad  =  all  —  0h  (42) 

The  maximal  invariant  set  C  within  Ee  is  where  equa¬ 
tion  (41)  and  (42)  are  satisfied. 

Notice  that  although  we  have  introduced  two  unknown 
variables  a  and  (3,  we  can  still  solve  for  (Zl5  Z2,  Ai,  A2) 
in  terms  of  {Id,  Ad,  SI,  8 A)  because  we  have  two  “extra” 
equations 

h-Ai=  0  l2  ■  A2  =  0  (43) 

Let  Ci  denote  the  set  where  a  =  0  and  0  =  0.  Let  C2 
denote  the  set  where  a  =  0  and  (3  ^  0.  Let  C3  denote 
the  set  where  a^0  and  0  =  0.  Let  C4  denote  the  set 
where  a/0  and  /3  0.  Then 

C  =  Ci  U  C2  U  C3  U  C4  (44) 

Among  all  the  possible  invariant  sets  we  have  calcu¬ 
lated,  only  Ci  is  the  one  we  want  the  system  to  ap¬ 
proach.  Thus,  we  shall  pick  suitable  values  for  SI,  SA, 

Id,  Ad  and  initial  conditions  so  that  Ci  is  the  only  pos¬ 
sible  invariant  set  within  H.  The  following  proposition 
gives  a  set  of  sufficient  conditions  to  achieve  this  goal. 

Proposition  3.5  Let  SI,  SA,  Id, Ad  satisfies  all  the  as¬ 
sumptions  in  proposition  (3.2)  and 

81  ■  5  A  =  0 


SI  ■  (Id  -SI)  =  0 
8 A  ■  (Aa  -  SA)  =  0 

(ld  -  SI)  ■  (Ad  -  SA)  =  0  (45) 

and  one  of  the  following  inequalities 

Id  7^  SI  —  aSA 

Ad  =£  SA  —  a  81  (46) 

where  a  =  ±\/2. 

Select  the  initial  condition  (qo,Po)  s.t 

V (qotPo)  <  c  (47) 

where  c  is  defined  in  proposition  (3.2).  The  system  will 
be  controlled  to  the  set  C\  where 

l\  —  I2  —  SI  =  0 
Ai  —  A2  —  SA  =  0 

h  =  Id 

Ai  =  Ad  (48) 

are  satisfied. 

Proof:  We  will  show  that  the  maximal  invariant  sub¬ 
set  of  the  compact  set  determined  by  the  given  initial 
conditions  contains  only  the  set  C\.  Other  invariant 
sets  will  vanish  or  will  be  impossible  to  reach. 

On  C2,  the  following  equations  are  satisfied 

h-l2  =  SI 

A!-A2  =  SA  (49) 


1 1  —  a.A\  —  Id 

Ai  —  ali  =  Ad  (50) 

from  (50),  since  ld  ■  Ad  =  0  and  l\  ■  A\  =  0  we  have 

a(\\hf  +  \\Ai  ||2)  =0  (51) 

Because  a  7^  0,  this  implies  that  l\  and  Ai  vanish. 
Thus  equation  (50)  will  not  be  satisfied.  C2  vanishes. 

On  C3,  the  following  equations  are  satisfied 

l\  —  I2  =  SI  +  /3A2 
A\  —  A2  =  SA  +  (3l2  (52) 


h  +  (3A2  —  Id 

A\  +  (3l2  =  Ad  (53) 

From  (52)  we  can  solve  for  1 1  and  A\  and  put  them  in 
(53).  We  get 

l2  +  ‘2,j3A2  =  Id  —  SI 

A2  +  2(3l2  =  Ad  -  SA  (54) 

Since  (ld  —  SI)  ■  (Ad  —  SA)  =  0  and  l2  ■  A2  =  0  we  have 

2/3(||  l2  ||2  +  ||  A2  ||2)  =  0  (55) 


Because  /3  ^  0,  this  implies  that  l2  and  A 2  vanish. 
Hence  C3  vanishes. 

On  C4,  the  following  equations  are  satisfied 

l\  —  l2  =  SI  +  (3A2 
A\  —  A2  =  SA  +  pi2  (56) 

li  —  aA\  +  [3A2  =  ld 
A\  —  ali  +  /3l2  =  Ad  (57) 

Replace  (3A2  and  (3l2  of  equation  (57),  we  have 

2/i  —  l2  —  aAi  =  Id  +  SI 
2Ai  —  A2  —  ali  =  Ad  +  SA  (58) 

From  (56)  we  can  solve  for  1 1  and  Ai  and  put  them  in 
(58).  We  get 

(1  —  a/3)l2  +  (2/3  —  a)A2  =  Id  —  SI  +  aSA 
(2/3  —  a)l2  +  (1  —  a(3)A2  =  Ad  —  5 A  +  aSl  (59) 

According  to  (45),  we  have 

(ld  —  SI  +  aSA)  ■  (Ad  —  SA  +  aSl)  =  0 

so 

(1  —  a(3)(2(3  —  or) (||  l2  ||2  +  ||  A2  ||2)  =  0 

Because  ||  l2  ||  and  ||  A2  ||  can  not  vanish,  there  are  three 
possibilities: 

(1)  cc  =  2/3  but  a/3  7^  1 
From  equation  (59)  we  have 

h  =  T^(ld-Sl  +  2pSA) 

^2  =  T-^(Ad-5A  + 2(351)  (60) 

Solve  for  l\,A\  from  equation  (56),  we  have 

h  =  ^  _  2/32  ^ d  +  fit^d  +  SA)) 

Ai  =  ^  _  2/32  (Ad  +  (3(ld  +  SI))  (61) 

Now,  since  li  ■  Ai  =  0  we  have 

(Id  +  (3(Ad  +  SA))  ■  (Ad  +  (3(ld  +  51))  =  0 
which  implies 

(3(\\ldf  +  \\Adf  +  \\Slf  +  \\SAf)  =  0 
This  result  is  impossible. 

(2)  a/3  =  1  but  a  7^  2/3 
From  equation  (59)  we  have 

l2  =  77-7 - (Ad  —  SA  +  aSl) 

2/3  —  a 
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Ai  —  — - (ld  —  51  +  a8A)  (62) 

2/3  —  a 

Solve  for  l\,A\  from  equation  (56),  we  have 

h  =  -zt. - (Ad  +  (3(ld  +  SI)) 

Ip  —  a 

Ai  =  -zt. - (Id  +  P(Ad  +  5A))  (63) 

2/3  —  a 

Now,  since  l\  ■  A\  =  0  we  have 

(Ad  +  (3(ld  +  SI))  ■  (Id  +  (3(Ad  +  <L4))  =  0 

which  implies 

0(\\  Id  ||2  +  ||  Ad  ||2  +  ||  SI  ||2  +  ||  SA  ||2)  =  0 
This  result  is  also  impossible. 

(3)  a/3  =  1  and  a  =  2/3. 

This  implies  that  a  =  ±\/2.  Equation  (59)  will  be 

Id  —  51  +  a5A  =  0 
Ad  —  5  A  +  aSl  =  0  (64) 

This  is  impossible  because  it  violated  the  conditions  in 
(46). 

Hence  C4  vanishes.  ■ 

A  special  case  of  this  proposition  is  when  we  choose 
51  =  0  and  5 A  =  0.  This  is  the  case  where  two  satellites 
are  driven  to  the  same  orbit. 

If  we  let  the  Lyapunov  function  V(J(q,p ))  = 

^[||  h  —  Id  ||  +  ||  A\  —  Ad  ||  ],  we  can  control  a  single 
satellite  to  transfer  between  elliptic  orbits.  This  case 
is  first  analyzed  in  [1].  In  this  case,  the  invariant  set 
is  simply  the  set  where  l\  =  ld  and  A\  =  Ad  are  satis¬ 
fied.  Comparing  to  this  case,  the  invariant  set  of  two 
satellites  case  is  more  complicated. 

To  completely  set  up  the  formation  completely,  we  need 
to  know  the  final  positions  of  the  satellites  on  their 
orbits.  This  can  be  computed  off-line  by  simulation. 
Then  we  can  choose  proper  starting  separation  between 
two  satellites. 


4  Simulation  Results 

In  order  to  verify  the  control  algorithm,  we  wrote  a 
simulation  of  controlling  two  satellites  into  formation 
on  MATLAB.  We  take  the  Lyapunov  function  as 

V(J(q,p))  =  ±(b1\\h-l2-5l\\2+ 
52||Ai  —  A2  —  <L4||2  +  &i||^i  —  ld ||2  +  &2IIA1  —  Ad\\2)  (65) 


here,  we  put  factors  61  and  62  into  the  Lyapunov  func¬ 
tion  so  that  the  control  is  weighted.  Our  control  is 

ui  =  —  £[&i(^i  —  h  ~  51  +  li  —  ld)'XQi 

+£1  x  62  (Ai  —  A2  —  5  A  +  A\  —  Ad) 

+(62(^1  —  A-2  —  <5A)xpi)xgi] 

U2  =  ?7[£>i  (Zi  —  I2  —  5l)xq2  +  ?2x^2(Ai  —  A2  —  5  A) 
+  (MAi  —  A2  —  5A)xp2)xq2]  (66) 

where  £  ,  77  are  positive  numbers  which  can  be  adjusted 
for  numerical  performance. 

In  practical  applications  this  control  law  needs  to  be 
discretized.  Here  we  use  a  simple  technique  to  ob¬ 
tain  a  discrete  control  law  from  (66).  Assume  that 
the  thruster  of  a  satellite  can  be  fired  towards  any  di¬ 
rection.  Let  T  be  the  time  interval  between  two  firings 
of  the  thruster  pulses.  Let  to  denote  the  starting  time 
of  the  controller.  Denote  Ui  the  discretized  control.  At 
time  to  +  nT ,  we  have: 

Ui(n)  =  Ui(t0  +  nT)  ■  T  (67) 

That  is,  we  assume  that  during  the  time  interval  T  the 
control  Ui  is  constant.  Of  course,  if  T  is  too  large  the 
algorithm  may  fail  to  converge. 


Figure  1:  Two  Satellites  Formation 

Figure  1  shows  when  51  =  0  and  <L4  =  0.  By  applying 
the  discretized  control  law  we  can  drive  two  satellites 
Satl  and  Sat2  onto  the  same  orbit.  We  take  the  unit 
length  to  be  the  radius  of  the  earth,  the  unit  time  to  be 
one  minute  and  the  unit  mass  to  be  1000kg.  The  ini¬ 
tial  conditions  of  the  two  satellites  are  given  by  specify¬ 
ing  their  six  orbital  elements  (a,  e,  i.  u>,  H,  r)  .  We  have 
(3,  0.3,  0,  tt/2,  0, 0)  for  Satl  and  (4, 0.2,  tt/4,  tt/4,  tt/3,  0) 
for  Sat2.  The  destination  orbit  is  a  circular  orbit  given 
as  (3,  0, 7r/6,7r/2, 0,  0).  Hence,  ld  and  Ad  are  deter¬ 
mined.  During  the  simulation  process,  we  noticed  that 
by  choosing  62  =  IOOO61  and  £  =  77  =  0.01  we  can  get 
a  decent  result. 


Figure  2  shows  the  change  of  ||  l\  —  h  ||,  100  ||  A\  —  A2  || 
and  the  Lyapunov  function  V  with  respect  to  time  dur¬ 
ing  the  whole  process.  As  we  can  see,  the  Lyapunov 
function  is  being  reduced  during  the  whole  process. 


xicr1 


Figure  2:  Above:  ||  h  —  h  ||  and  ||  Ai  —  A 2  ||  as  a  function 
of  time;  Below:  the  Lyapunov  function  as  a 
function  of  time 
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Figure  3:  The  final  separation  between  two  Satellites  as 
a  function  of  starting  separation 

The  final  relative  position  of  the  two  satellites  depends 
on  the  starting  relative  position  of  the  two  satellites. 
In  the  following  experiments,  two  LEO  satellites  are 
first  placed  on  the  same  circular  orbit  specified  as 
[1.1, 0,  0, 7r/2, 0,  0]  with  certain  amount  of  initial  sep¬ 
aration.  Then  they  are  transfered  to  a  final  orbit  of 
[1.3,  0, 0,  7t/2,  0,  0].  Figure  3  shows  the  final  separation 
versus  the  initial  separation  between  the  two  satellites. 
It  can  be  seen  that  the  separation  is  increased  since  we 
are  transferring  to  a  higher  obit. 


